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Dynamical twisted mass fermions 



1. Introduction 



Within the framework of Symanzik's improvement programme [Q, ||] the emphasis in lattice 
QCD has been on the construction and testing of fermion actions to avoid 0{a) lattice spacing 
artefacts. Since the lattice artefacts induced by the gauge action start with 0{a^), the need of also 
exploring different gauge actions has not been considered as equally urgent. However, although 
this point of view is certainly valid in the asymptotic regime of small lattice spacing a, the situation 
is different when values of the lattice spacing, say, a > 0. 1 fm are considered. In this region higher 
order lattice artefacts my play an important role. 

In the case of Wilson fermions the explicit breaking of chiral symmetry leads to the appear- 
ance of a strong first order phase transition that severely affects the numerical simulations [^, |8|]. 
As we will demonstrate in this contribution, the choice of the gauge action can have a strong effect 
on the strength of this phase transition and therefore it becomes important to consider both, the 
fermion and the gauge action to find a suitable lattice QCD action for simulations of dynamical 
quarks. 

The twisted mass fermions approach [^] provides, among other advantages, the ideal frame- 
work for the investigation of the zero-temperature phase diagram of lattice QCD with Wilson 
fermions, see refs. [10, 11, 12] for reviews on twisted mass fermions in present and past con- 
ferences. 

Our present understanding can be summarized as follows, see also fig. [l[ For values of the 
lattice spacing much coarser than a = 0.15 fm an unphysical phase appears for small quark masses, 
the Aoki phase |13, 14, 15]. The transition from the standard lattice QCD phase to the Aoki phase 
is of second order: here the charged pions become massless. The Aoki phase only extends in the 
untwisted quark mass direction being absent for non vanishing twisted quark mass pL. For smaller 
values of the lattice spacing a first order phase transition appears |8|, [l^, |l^, this time extending 
in the twisted mass direction pL. The phase transition occurs at vanishing (untwisted) quark mass 
and separates the phases with opposite signs of the quark mass. This first order phase transition 
is reminiscent of the continuum phase transition when the quark mass is changed from positive 
to negative values with the corresponding jump of the scalar condensate as the order parameter of 
spontaneous chiral symmetry breaking. The generic phase structure of lattice QCD was discussed 
in refs. [0,|,|l6|]. 



The appearance of the first order phase transition has serious consequences, since in such a 
scenario the pion mass cannot be made arbitrarily small but assumes a minimal value, m™^, 
which may be about 500 MeV and hence it becomes impossible to work close to the physical 
value of the pion mass. It therefore becomes important to understand the phase structure of lattice 
QCD as a pre-requisite before starting large scale simulations. Our collaboration has performed a 
detailed study of the phase diagram of lattice QCD. We used Wilson fermions with and without 
twisted mass pai^ameter and studied the phase diagram as a function of the lattice spacing for the 
Wilson plaquette, the DBW2 and the tree-level Symanzik improved gauge actions. As we will 
show in this contribution, the strength of the first order phase transition, namely the size of m^'", 
depends strongly on the choice of the gauge action when comparable physical situations are tested. 
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Figure 1: Current knowledge of the lattice QCD zero-temperature phase diagram with Wilson fermions as 
a function of the inverse gauge coupling j5 °^ 1/^^' 'he hopping parameter K and the twisted mass parameter 



Since lattice chiral perturbation theory (;t;PT) |^ [18 , 19, 2^ 21 , 22] predicts a weakening of the 
first order phase transition toward the continuum limit, it is interesting to check this prediction and, 
in particular, to investigate how fast the transition weakens when the continuum limit is approached. 
The answer to the latter question will naturally depend on the choice of the actions that are used 
for the gauge and the fermion fields. Moreover, the predictions of ;^PT for e.g. the quark mass 
dependence of the pion mass and the pseudoscalar decay constant at such a phase transition can be 
directly confronted to results from numerical simulations. 

This contribution is organized as follows. In section 2 we will give the definitions of the 
actions we will use. In sections 3, 4 and 5 we will give our results for the Wilson plaquette, the 
DBW2 and the tree-level Symanzik improved actions, respectively. In section 6 we will confront 
our lattice data with lattice ;tPT. In section 7 we will provide first results from simulations with 
Nf = 2+1 + 1 flavors of quarks. Finally, we will conclude in section 8. 



2. Lattice action 

The lattice action for a doublet of degenerate twisted mass Wilson fermions (in the so-called 
"twisted basis") reads 




(2.1) 



with = + 4r = 1/2k, r the Wilson-pai^ameter, set in our simulations to r = 1, aniQ the bare 
"untwisted" quark mass in lattice units (k" is the conventional hopping parameter) and jj. the twisted 
quark mass; we also define Ux-^ = U^ n „ and = —Yn- 
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For the gauge sector we consider the one-parameter family of actions including planar rectan- 
gular (1x2) Wilson loops (U^^^): 

Ss = P'L(co I \i-\rcU^^A+c: £ |l-^Ret/^x,4V (2.2) 

X \ ^<v;^,v=l I -* J p^v;^,v=l I J J 

with the normalization condition cq = I — 8ci. We consider in this contribution the three cases: i.) 



Wilson plaquette gauge action, ci = 0, ii.) DBW2 gauge action [|23[], ci = —1.4088, iii.) tree-level 
Symanzik improved gauge action (tlSym) [El], ci = — 1/12. 



3. Wilson plaquette gauge action 

The first action we investigate here is the Wilson plaquette gauge action. We studied the lattice 
spacing dependence of several quantities [|l^] keeping the twisted mass and the lattice size roughly 
fixed to roH ^ 0.03 (ro being the Sommer parameter [25], we assume ro = 0.5 fm throughout this 
contribution) and L f» 2 fm. We varied j8 = in the range j8 = 5.1 — 5.3. For each value 
of the hopping parameter we performed a hot and a cold start in order to check for co-existing 
values of physical observables which we have chosen as the pion mass niji and the rho meson mass 
jiip, extracted from the usual correlation functions. Another quantity is the amplitude f^^, which 
reduces to the pion decay constant /;;: for /X = 0, defined as 



fy = — P(0 71 = — gj, , where vap - ' ' 



(3.1) 



(0|P(0)|7r) 

is extracted from the asymptotic behavior of the ratio of the axialvector-pseudoscalar correlator 



Cap{xq) to the pseudoscalar-pseudoscalar correlator Cpp{xq) (cf. ref [ |l7[ ] and references therein 
for details): '^'^Jl}^°J^ = ?"Aptanh [mji{T /2 — xo)] . Moreover we measure the untwisted PCAC quark 



mass m 



PCAC. 



m 



PCAC 
X 



2(0|P(0)j7r) 2g^ 



(3.2) 



Note that the physical decay constant (PCAC quark mass) is obtained by a combination of f^^ and 
the twist angle CO {m^y^^'^ and the twisted mass /i), see also section 4. 
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Figure 2: 

mass (right). 
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Lattice spacing dependence of the plaquette expectation value (left) and the PCAC quark 
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The plaquette expectation value shows the typical behavior at a first order phase transition 
with meta-stable branches, illustrated in fig. ^ (left). For the PCAC quark mass the two branches 
correspond to positive and negative quark mass, respectively (see fig. § (right)). In both cases the 
gap decreases from /3 = 5.1 to 5.3 and the meta-stability region in 1 /2k shrinks. For the pion mass 
shown in fig. || (left) this fact has the important consequence that only a minimal pion mass 500 
MeV can be reached. If one attempts to lower the mass further by tuning jj. or k:, a jump to the 
other phase occurs. In addition, the scale parameter ro/a is quite different in the two phases and its 
mass dependence is non-negligible, as illustrated in fig. ^ (right). 
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Figure 3: Pion mass vs. PCAC quark mass (left) and mass dependence of ro/a (right) at /3 = 5.3. "Low" 
and "high" refer to the low plaquette (positive quark mass) and the high plaquette (negative quark mass) 
phase, respectively. 



In the light of these observations the interesting question arises at which j3 -value or lattice 
spacing we can safely simulate pions with masses as light as 300 MeV in order to make contact 
with ;tPT. With the present data it is very difficult to make an extrapolation, since there is a large 
ambiguity in the determination of the minimal pion mass, the lever arm is very short due to the small 
range of j3=5. 1-5.3 covered and furthermore due to the large difference of the scale parameter r^/a 
in the two phases. Nevertheless, a qualitative estimate would be a lattice spacing of 0.07-0. 1 fm, 
where one should be able to reach small pion masses without being affected by the first order phase 
transition. This makes even L = 2 fm simulations, e.g. for a detailed scaling study, very demanding. 
Therefore alternative gauge actions that lead to a reduced strength of the first order phase transition 
are investigated in the next sections. 



3.1 Scaling behavior 

Here we study the scaling behavior with our current data for the Wilson plaquette action, also 
including the data from the DBW2 gauge action, see sec. ^. To this end we express the physical 
quantities in dimensionless variables. We first define a reference quai^k mass by {rotn^Y = 1.5. At 
this reference quark mass, we then determine* reference values for the quantities we are interested 
in, i.e. ro/a, m^^^^, rriji, trip and f^^. This allows to define dimensionless ratios such as for 

the quark mass, (7 — pcyi(; , — or other observables, Rq — ^ P . The so determined observables, 

^^^X I ref ' l ef 

measured in units of their reference values, are universal and can be compared at different j8 -values 
and for different actions. 

*We consider only the positive quark mass phase which corresponds to standard lattice QCD. 
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Figure 4: Scaling ratios R^ps in the quenched approximation (upper panel) and for dynamical 
fermions (lower panel). The lower panel also includes results from the DBW2 action at j3 = 0.67 and 
j3 = 0.74 (see also sec. 4). 



In fig. (lower panel) we show as an example the scaling ratio for our dynamical fermion data 
of the Wilson plaquette and the DBW2 gauge actions. We observe no scaling violations for the 
data from the Wilson plaquette gauge action at j3 = 5.1 — 5.3 within the (large) statistical eiTors. 
Moreover, we see agreement with data at j8 = 5.6 for pure Wilson fermions without twisted mass, 
produced for the algorithmic study of an HMC variant [26, 27]. For the data from the DBW2 gauge 
action at smaller masses we observe instead a slight difference in the scaling behavior of the n = 
and iU 7^ results. This might indicate that such a difference could be detected for the Wilson 
plaquette gauge action only at higher beta values, where smaller quark masses can be reached. 
Similar plots were also obtained for the scaling ratios /?„,2 and niji/mp. For comparison, we also 
show in fig. ^ (upper panel) an example of the ratio R^ps for quenched data at full twist, K = Kc [p^]. 
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Here we observe scaling violations for our coarsest lattice at j3 = 5.7 and a difference in the scaling 
behavior at small masses for different discretizations of twisted mass fermions, realized by using 
two definitions of the critical mass {Kc'"" and K^^^^ , see ref. [^]). This might provide a warning 
that at small values of the quark mass the lattice artefacts could be significant and that scaling 
violations might show up. 

4. DBW2 action: analysis of twist angle and physical quantities 

In this section we analyze several interesting quantities in the context of simulations of Nf = 2 
flavors of twisted mass fermions with the DBW2 gauge action. This particular setup was studied 
in [16] on a 12^ -24 lattice at j3 = 0.67 with spatial extension L ~ 2 fm, lattice spacing a ~ 0.18 fm, 
ro/x ~ 0.03: the minimal pion mass found here was lUji ~ 360 MeV. Signs of metastabilities ai^e 
hardly detectable in this situation, which is a clear improvement compared to the Wilson plaquette 
gauge action, cf. sec. 3. In this section we also consider data from 16^ • 32 lattices at a finer value 
of the lattice spacing, j3 = 0.74, keeping approximately the same volume and value of the twisted 
quark mass as in the aforementioned simulations. 

4.1 Twist angle 

The twist angle ft) defines the chiral rotation relating twisted mass QCD to ordinary QCD. In 
the case of the chiral currents the rotation reads (considering only charged currents, a= 1,2): 

V"^ = ZyV"^ COS (o + Sab ZaA''^^ sin (0 , (4. 1) 

A"^ = ZaA"^ cos (0 + eahZyV^^ sin co (4.2) 

where the hatted currents on the l.h.s. denote the chiral currents of QCD (physical currents), while 
the cun^ents on the r.h.s. are the corresponding bilinears of the quai^k-field in the twisted (x-) basis. 
Note that the renormalization constants of these bilineai^s, Zy and Za, ai^e involved. For a given 
choice of the lattice parameters, the twist angle (O is determined by requiring parity conservation 
for matrix elements of the physical currents [||, 16 1. Since unknown renormalization constants are 



involved, two conditions are required, our choice being: 

£(y+p-)=o, (4.3) 

I(Ain7)=0. (4.4) 



The solution of eqs. (4.3) and (4.4) with eqs. (4.1) and (4.2) gives a direct determination of the 



twist angle (O and of the ratio Z^/Zy from lattice data, see [ |16| ] for details. Eq. ( t4.3| ) implies e.g. 



cot ftv ^ ^ cot co = -i . (4.5) 

In particular at full twist where (O = n/2 the condition reads ^(^^^r) = or cotftv = 0. The 

X 

full twist situation can be also obtained by requiring the vanishing of the PCAC quark mass m^J^^^. 
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One can easily see that using the PCAC quark mass as a criterion for full twist just amounts to a 
different definition of the twist angle: the current is replaced by its divergence in (p3|). So the two 
criteria ai^e equivalent and differ (at most) by 0{a) effects. The numerical equivalence of the two 
criteria for full twist is illustrated by fig. ^ 




1/2K 1/2k 



Figure 5: Determination of the critical hopping parameter Kc corresponding to full twist by parity -restoration 
and by extrapolating to zero the untwisted PCAC quark mass m^^^^ . The small discrepancy observed at 
/3 = 0.67 (left panel) between extrapolations from positive and negative quark masses is probably a residual 
effect of the first order phase transition (crossover). For j3 = 0.74 (right panel), extrapolations from both 
sides give consistent results. 



4.2 Physical quark mass and pion decay constant 

The knowledge of the twist angle is necessary for the determination of physical quantities like 
the quark mass and the pion decay constant. The physical PCAC quark mass m^'"^^ can be obtained 
from the Ward identity for the physical axialvector current: 

(V;A+/^7) = 2m™(P+P7). (4.6) 



From eqs. ( |4.1[ ) and ( |4.2| ) follows 

^"ii = ^xii cot (0 + Sab V^^ (sin (o)-^ (4.7) 

where the conserved vector cuiTcnt of the ;^-fields has been now considered, for which Zy = I. 
Inserting the above result in the Ward identity ( pT^ ) and using parity restoration for the physical 
currents one obtains: 



m^^^^ = -/(2sin«)-'— ii^^^ =ZA(cosw)-'m^^^^ . (4.8) 

{Px Py ) 

Analogously, for the physical pion decay constant fj^ we use 

f^ = m^'{0\A+{0)\7i+) = -/Ksin«)-i(0|yo+(0)|7r+) . (4.9) 
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Table 1: Chiral extrapolation (m^ = 0) of the Sommer scale parameter ro and 
pion decay constant f^; the scale independent combination fjir^ is also reported. Only 
data with positive twisted quark masses have been used for the extrapolations, with the 
exception of the point near full twist at jU = 7.5 • 10^^. 







ro/a 


a (fm) 


a fit 


fnrQ 


0.67 
0.74 


1.0-10-2 
7.5 - 10-3 


2.680(68) 
4.11(13) 


0.1866(72) 
0.1216(39) 


0.1171(59) 
0.0726(25) 


0.320(10) 
0.309(15) 



In fig. ^ the pion decay constant is plotted as a function of the quark mass (the simulation points 
for negative quark masses are not included). The figures also show the determination of fj^ by 
the axial vector current A"^: a formula similar to eq. (4.9) applies in this case, where however 
the factor 1 / sin ft) is replaced by 1 / cos ft). In the interesting region near full twist this factor 
introduces the large fluctuations in the estimate of fji observed in the data. Moreover in the case 
of the decay constant has not yet the right normahzation of the continuum: a Za factor is 
missing. In the case of the conserved vector current on the contrary the lattice determination 
of has automatically the correct normalization [29, 30]. If we exclude the lightest point at 
j3 = 0.67, which is likely to be under the influence of residual metastability effects, fj^ seems to 
be characterized by a linear dependence upon the quark mass. On the basis of this observation we 
try a simple linear extrapolation to the chiral limit rrf^^^ = 0, see table |l| for the numerical results. 
(Of course deviations from this linear behavior could be present for lighter quark masses, where 
chiral logarithms might play a role.) We see that the two j3-values give compatible values for fnro, 
suggesting that there are no large scaling violations for this action. These values, obtained with 
Nf = 2, are also near to the physical value /^ro = 0.33. 
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Table 2: The renormalization constants Zy, at full twist (extrapolating from posi- 
tive quark masses) with comparison to 1-loop perturbative estimates (PT) and tadpole- 
improved perturbative estimates (TI)[^. 



/3 








Op. 


Z 


Z(PT) 


Z(TI) 


0.67 


1.0 


•10 


-2 


V 


0.5650(11) 


0.6089 


0.6531 


0.74 


7.5 


•10 


-3 


V 


0.6217(23) 


0.6459 


0.6892 


0.67 


1.0 


•10 


-2 


A 


0.952(30) 


0.7219 


0.7176 


0.74 


7.5 


•10 


-3 


A 


0.944(74) 


0.7482 


0.7735 



4.3 Renormalization constants 

The renormalization constant of the vector current Zy can be determined on the basis of the 



non-renormalization property of the conserved current Vx^ [pip. In the case of twisted mass QCD 
we use 

Note that in twisted mass QCD the time-component of the vector current couples the vacuum to the 
pseudoscalar meson: in the most interesting region near full twist this coupling is (in our twisted 
basis) maximal. In lattice QCD with ordinary Wilson fermions the analogous procedure has to rely 
on the noisier matrix element of the vector particle. Eq. ( |4.10| ) represents a mass dependent pre- 
scription. We obtain a mass independent determination of Zy by extrapolating Zy" to full twist. In 
this situation the theory is 0{a) improved and Zy" has only 0{a^) lattice deviations. (In particular 
residual 0{{pLaY) terms can be eliminated by a — extrapolation.) The results are reported 
in table § the renormalization of the axialvector cuiTcnt Za is also reported: this is obtained by 



combining the direct determination of Zy ( [4.10| ) and the determination of Z^/Zy from the parity 



restoration (see above) after extrapolation to full twist. Observe that the conditions (43), ( p^ 
hold in general up to 0{a) parity violations: imposing (as we do) exact parity restoration for the 
considered m.e. leads to a mass dependent determination of the ratio Z^/Zy. 

5. Tree-level Symanzik improved gauge action 

In this section we will now concentrate on the tree-level Symanzik improved (tlSym) gauge 
action for which the coefficient c\ = — 1/12. The reasons for the choice of this action are manifold 
- firstly, the action is designed to show a good and smooth behavior in the perturbative regime; 
secondly, in many quenched scaling studies the action has proved to behave well also in the non- 
perturbative regime; and thirdly, the additional coefficient c\ is rather small in absolute value as 
compared to other improved actions like the DBW2. We have simulated several values of k close 
to Kc, i.e. close to maximal twist, at three different values of j3 using a HMC algorithm variant with 



mass preconditioning and multiple time scale integration [g^, p7| ]. 

In order to check for a possible phase transition and corresponding metastabilities we measure 
the average plaquette value as a function of the hopping parameter K on runs with hot (disordered) 
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and cold (ordered) starting configurations. Since the metastability, if any, will show up around 
we concentrate our hot and cold runs on k:- values closest to only. Fig. (left) shows the results 
for the average plaquette value for our coarsest lattice at j3 = 3.65, 12^ • 24 and roH = 0.038. Here 
and in the rest of the section the red circles and blue squares represent the results from the hot and 
cold starts, respectively. From fig. ^ (left) it is clear that there seems to be no metastability as far as 




K MC time 



Figure 7: Left panel: Average plaquette value vs. k at p = 3.65, ro/i = 0.038 on a 12^ • 24 lattice from 
hot (red circles) and cold starts (blue squares). Right panel: Plaquette MC time history for two runs at 
j3 = 3.65, ro/X = 0.038, k — 0. 17024 on a 12^ • 24 lattice starting from hot (red line coming from below) and 
cold configuration (blue line coming form above). 

the plaquette value is concerned, but the existence of a phase transition at this parameter set can of 
course not be excluded. Indeed, the rise of the plaquette value around JQ- ~ 0. 17025 points toward 
the fact that we ai^e at least in the vicinity of a phase transition and we conclude that a lower value 
of /I would hardly be possible to simulate in practice. The vicinity of the phase transition is also 
reflected in the MC time histories of the plaquette value. We observe a rather strong critical slowing 
down leading both to a very slow thermalization and large fluctuations of the plaquette value over 
several hundreds of trajectories. It is evident from the MC time history shown in fig. |7| (right) that 
any statement about the plaquette value itself as well as any eiTor estimate are very difficult to make 
and should hence be taken with care. The figure also illustrates the fact that simulations close to 
or at Kc (i.e. at maximal twist angle) are rather difficult, at least at such large values of the lattice 
spacing (a ~ 0.13 fm) and values of jJ. as chosen here. 

Let us now move to our simulations at j3 = 3.75. Fig. |8] (left) shows the average plaquette 
value vs. K" at j3 = 3.75, ro/x = 0.020 on a 12^ -24 lattice from hot (red circles) and cold starts 
(blue squares). Compared to fig. we are now at finer lattice spacing {ro/a = 3.76(9) at j8 = 3.65 
vs. ro/a = 4.1(2) at j8 = 3.75), but on the other hand also our value of the twisted mass is halved, 
so we do not necessarily expect the situation to be better. Indeed, we find that at the value of 
K = 0.1660 that is closest to JC^ the average plaquette value from the hot and the cold start do 
not match each other, which is pointing to the existence of a phase transition and a coiTcsponding 
metastability at this set of parameters. On the other hand, considering the large fluctuations of 
the plaquette similar to the ones in fig. ^, the runs here are presumably not long enough to allow 
reliable estimates and the metastability might disappear as the MC time history becomes longer. 
Nevertheless, the strong rise of the plaquette around jq. points toward the fact that we are at least 
close to a transition, even though the rise appears to be much weaker than at jS =3.65. 
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Figure 8: Left panel: Average plaquette value vs. K at /3 = 3.75, ro/i = 0.020 on a 12^ • 24 lattice from hot 
(red circles) and cold starts (blue squares). Right panel: Average plaquette value vs. K" at j3 = 3.90, rg/i = 
0.041 on a 16^ • 32 lattice from hot (red circles) and cold starts (blue squares). 



The situation at j3 = 3.90 is illustrated in fig. |8| (right) where rofX = 0.041 corresponds again 
to our larger value of /i. Here we observe a very smooth dependence of the plaquette value on K 
and there seems to be no trace of a cross-over or even a nearby phase transition. (We attribute the 
fact that the plaquette value at k: = 0.1610 ~ Kc from the hot and cold start do not match exactly to 
the shortness of our runs which prevents an accurate determination of the average plaquette value.) 



Next we consider the untwisted PCAC quai^k mass m^'-^'^, cf. eq. (3.2), as a function of K for 



j8 = 3.75 and 3.90 in fig. For /3 = 3.75 (left) we find that the mismatch between the plaquette 




Figure 9: Untwisted PCAC quark mass m^-^ vs. K from hot (red circles) and cold starts (blue squares). At 
j8 = 3.75, roM =0.020 on a 123-24 (left plot) and j3 = 3.90, rojU = 0.041 on a 16^-32 lattice (right plot). 

value from the hot and the cold start at k: = 0.1660 ~ is also reflected in the PCAC quark mass, 
although we expect that this mismatch will eventually go away for long enough runs as discussed 
before. Nevertheless, even then it seems unlikely that extrapolations to m^'-^'' = from K"-values 
above and below jq. will coincide. This is in contrast to the situation at j8 = 3.90 (right plot) where 
we observe a smooth dependence of m^^^^ on K throughout a wide range of KT-values. Indeed, 
extrapolations to m^'^'^'^ = from K"-values above and below Kc nicely coincide while the slopes 
are different (in consistency with expectations from ;tPT [18, |l^, 20, 21, ^]). 

Let us now turn to the pion mass squared as a function of the PCAC quark mass tti^j^^^ 
which we plot in fig. pi for j8 = 3.75, ro/x = 0.020, 12^ • 24 on the left and for /3 = 3.90, ro/x = 
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0.041, 16^ • 32 on the right. At j8 = 3.75 the pion mass that is realized at roH = 0.020 is nijc 
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Figure 10: Pion mass squared vs. PCAC quark mass mx^'~^'^ from hot (red circles) and cold starts 
(blue squares). At j3 = 3.75, rojU = 0.020 on a 12^ ■ 24 (left plot) and j3 = 3.90, rojU = 0.041 on a 16^ ■ 32 
lattice (right plot). 

353(22) MeV for the hot start and = 382(33) MeV for the cold start, so ~ 370 MeV seems 
to be a good estimate for the minimal pion mass that can be reached at this parameter set. For 
additional simulations at the same j3 but on a larger volume 16^ • 32 at k: ~ jq. and tq/x = 0.039 we 
obtain = 407(16) MeV. Finally, at j8 = 3.90 the pion mass we reach at roji = 0.041 is ttiti = 
453(32) MeV, while an additional simulation at K" ~ and r^pL = 0.020 gives nijc = 279(37) MeV 



6. Chiralfits 



The extension of ;^PT to the case of adding a twisted mass term was considered in refs. [33 



18, 1^ 20, 21, 22]. For the low-energy constants in next-to-leading order [34, Q we use the no- 
tation L54 = 2L4 +L5, L86 = 2U+U, = 2W4 + W5, ^86 =2W(,+W^,W = \{Wif, - 2L86), 
W = 5(VK54 — L54) where the L, ai^e the Gasser-Leutwyler coefficients and the Wi are the corre- 
sponding coefficients of the 0{mga) coiTcctions. We consider the dependence of lattice quantities 
on the untwisted PCAC quark mass m^^^'-- , so we define XpcAC ~ ^i^^r"^ where we explicitly 
express that m^^^ is renormalized in some prescription: m^^j^^ = ZAZp^m^^^ . We also define 

the combination containing the physical quark mass % = IBq^J {m^^^Y + M| with jiR = Zp^fj.. 
Defining p = 2Woa, where Wo is the low-energy scale for 0{a) breaking (analogous to Bq), Fq the 
pion decay constant in the chiral limit, and A the renormalization scale, ;tPT at next-to-leading 
order with 0{a) corrections gives: 

=^ + ^ij2^X''";^ + |2((-'-M + 2'-M)«' + 2("'-»')PXrac} ((5 1) 

t-'-16^^'4 + 4<^"^ + ^*''^' <«' 
H = >-n^«l"| + 4((-^« + 4.,a)x + (4«/-2^)p^|^). (6.3, 

We set A = AtiFq. 
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We performed a combined chiral fit of tlie quantities m| , and gjc for our data from tlie Wilson 
plaquette and the DBW2 gauge action (including three and two different j8 -values, respectively, see 
sec. 3 and 4) in order to obtain estimates for the combinations of coefficients Lse, ^^54, W, W and 
of the low-energy constants Bq and Fq. In fig. [ll] we show the result of this global chiral fit for the 
case of the pion mass, for the Wilson plaquette action (left) and the DBW2 action (right). In the 
former case three values of p are included. 



Pion Mass, p = 0.67, anj = 0.01 




Figure 11: Chiral fit of the pion mass; Wilson plaquette action (left) and DBW2 action (right). In the Wilson 
case a darker blue tone corresponds to a larger value of )3 . 

The general observation is that lattice ;^PT theory formulae reproduce rather well the lattice 
data. The 0{a) terms of the chiral expansion are in general not dominant in the fits, which one can 
infer from the fact that the associated coeffiecients fluctuate with large deviations around zero. This 
is probably due to the use of the PCAC quark mass which in past experiences has shown to have 
reduced lattice artefacts compared to the I /2k definition of the bare quark mass. In fact fig. § (right 
panel) shows the presence of lai^ge lattice artifacts. 

The use of different setups for the fits (different sets of data-points and actions, independent fit 
methodologies) allowed to test the stability of the results for the low-energy constants. Combining 
aU information, we get the following ranges (here we fix the lattice scale by using the chiral- 
extrapolated value of the Sommer scale ro): Fq = 70 - 100 MeV, Bo/Zp = 3-6 GeV, L54 = (0.8 - 
1.8) • 10^^, L86 = (0.5 — 1.0) • 10^^. As mentioned, the coefficients of the 0(a) terms could not be 
determined. 



7. Nf = 2 + 1 + 1 flavors of twisted mass fermions 

In this section we consider simulations of QCD in the twisted mass setup including the dy- 
namics of the strange (s) and charm (c) flavors. The u and d quarks, which are degenerate, are still 



described by the action ([2.1|). 



7.1 Split doublet 



The mass-splitting in the charm-strange (cs) sector is realized by adding an extra mass-term [ 36] 
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in the action for a doublet of degenerate twisted mass fermions: 



where we use the same notation as in eq. (2.1 ), in particular jj.^^ = aniQ + 4r = 1 /2Kcs- The new 
term with coefficient /ig (the minus sign is conventional) is responsible for the splitting. This 
construction preserves reality and, for jX^^ > jjLg, positivity of the fermionic determinant. The 
particle content of the theory is more transparent in a different basis when Xx ^^p{ij'^2}Xx^ 
Xx Xx^'>'-p{—i^'^2}, where the mass term in the untwisted direction is diagonal and the twist is 
in the Ti flavor direction: 

Sq,ma^, = Y^Xxi^K + ^sa^S + m^iaiXcsjXx ■ (7-2) 

X 

In analogy with the degenerate case the "average" quark mass of the split-doublet is given by (we 
neglect here for simplicity renormalization factors) ni„ = {m^' — wIq*.^)^ + jJ-cs^ and the individual 
masses by tric^s = ± Ms- At full twist: 

niu = nid = l^ud , rric = lies + l^s , m= l^cs - l^s ■ (7-3) 

7.2 Tuning to full twist 

Tuning to full twist the Nf = 2+l + l theory looks apparently more complicated because of 
the presence of two tunable variables and Mq' (or hopping parameters K"„(/ and K"„) with their 
respective critical values. This problem can be however circumvented in a strategy analogous to 
the Nf = 2 case. 

As we have seen previously in this contribution for the Nf = 2 theory, one possible definition 
of the critical quark mass mo„{gQ,iJ.) is given by the vanishing of the corresponding PC AC quark 
mass m^^'-- . Due to chirality breaking the latter gets shifted: 

m^^^'^ = mo - a^^ figQ,amo,an) , (7.4) 

with / a dimensionless function. On the basis of the symmetry of the action under parity x (/i ^ 
—jj.) one can show that the additive renormalization of the quark mass is even in jj., and analyticity 
in turn implies 

f{go,amo,an)=f{gQ,amo) + 0{lJ.'^a^) , (7.5) 

where f{go,amo) is the shift for ordinary Nf = 2 QCD without twisted mass term. So the twisted 
mass term in the action only produces an 0{a) effect on the quark mass (with go and nio held fixed): 

m^^^^ = mo-a-^f{go,amo)+0{a) . {1.6) 

The above argument can be easily generalized to the Nf = 2 + \ + \ theory. Here one has to make 
a distinction between the two sectors: 

"^/M = - a^^ fud{go,am'(f ,am^Q\aHud,aHcs,aHs) , (7-7) 



PCAC 



"^yM = mo - a ^ fcs{gQ, ami' ,am'^' ,alXcs,allud, alls) . (7.8) 
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The functions /„j and /^j are in this case even in /i„t/, /x„ and /ig^: similarly to Nf = 2, the 
associated terms in the action only affect the additive renormalization of the quark mass by 0{a) 
terms. So we write: 

'"S = '"0^ - a-'f{go,am'^,am^') + , (7.9) 
PCAC _ a-'f{go,am:^,amf) + (9(a) , (7.10) 



where on the r.h.s. we have now the mass-shifts for the theory without twist and mass-splitting 



(Nf = 2 + 2 QCD): here the distinction between the two sectors is immaterial. From eqs. (7.S), 
( [7.10| ) it follows immediately 

< = < = '"o ^ m™ = mjc^^ + 0(a). (7.11) 

The above result suggests to tune rn'o' to its critical value where tn^^^^^ = keeping wiq* = tuq^ (or 
^cs = ^ud)'- in this situation m^^'^f = 0{a). Observe that since the average quai^k mass in the cs 
sector is typically large, the 0(a) error is expected not to affect the full twist improvement in the 



sense of [37], while it is critical to have good tuning in the light quark sector. 



7.3 Numerical simulations 

Our strategy is to simulate the Nf = 2 + 1 + 1 theory at full twist. As previously explained for 
a given value of jS this can be obtained by fixing K^., = K„d = K and tuning the PCAC quark mass 
in the light quark sector to zero (or equivalently by tuning the twist angle to n/2, see sec. At 
the critical value K = Kcr, '''^^^ud) = by definition and ni^^^^ = 0(a), where the 0(a) error in the 
cs sector is not expected to spoil 0(a) improvement. 

For the gauge action we choose the tree-level Symanzik improved action, which shows for 
Nf = 2 an acceptable behavior with respect to the phase transition at small quark masses and in 
addition presents several more advantages in comparison with the other possible candidate DBW2, 
see also section 5. 

At present, we are simulating on two lattice-sizes: 12^ • 24 and 16^ • 32, our preliminary ob- 
jective being the reproduction of the physical situation of previous simulation points [|l^] for the 
Nf = 2 theory. This means a lattice extension L = 2 fm and ro/a = 3 and 4 (coiTcsponding to lattice 
spacings a ~ 0.17 fm and 0.12 fm, respectively), Ughtest pion mass ~ 400 MeV (ni„^ ~ nis/S) 
and 300 MeV (m„d ~ wi,/5). 

The s and c quark masses need to be tuned to their physical values. To this end we measure the 
kaon and D meson masses and monitor the dimensionless quantities roniK and romo which in nature 
assume the values roniK — 1.25 and romo ~ 4.7. Due to the Ti term in the action, see eq. ([7^), 
kaons and D mesons get mixed in the twisted theory. We solve the problem by diagonalizing the 
two-by-two matrix of the K-D correlators. 

The analysis of the spectrum in the cs mesonic sector (D^) is even more awkward, again due 
to flavor mixing. The Ti term in the action introduces indeed off-diagonal matrix elements in the 

^An additional symmetry in the cs sector is needed for the argument, namely Xx ^ exp{;^Ti};)^j;, Xx ~> 
;f^exp{— i^Ti} composed with fig — M6- 
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Table 3: Simulation parameters of the performed simulations. See the text for the 
explanations of the symbols. 



run 


P 


K 








Nconf 


(a) 


3.30 


0.1720 


0.01 


0.325 


0.275 


550 


(b) 


3.30 


0.1725 


0.01 


0.325 


0.275 


920 


(c) 


3.30 


0.1730 


0.01 


0.325 


0.275 


900 


(d) 


3.30 


0.1735 


0.01 


0.325 


0.275 


730 


ie) 


3.30 


0.1740 


0.01 


0.325 


0.275 


590 



Table 4: Measured quantities: the Sommer scale ro; p meson mass mp, pion mass 
m,t and PCAC quark mass am^^^^ (light quark sector); kaon, D and meson masses. 
For the Dg meson only connected diagrams have been taken into account. 



run 


ro/a 




anip 


ariiK 


amj) 






(a) 


2.23(7)(10) 


0.7102(38) 


0.8884(41) 


0.9551(21) 


1.456(10) 


1.6456(49) 


0.08795(57) 


(b) 


2.40(5)(5) 


0.6069(61) 


0.8075(38) 


0.9078(41) 


1.420(11) 


1.6385(47) 


0.06399(73) 


(c) 


3.05(3)(8) 


0.3463(49) 


0.633(13) 


0.7940(16) 


1.247(19) 


1.5684(59) 


0.01589(82) 


(d) 


4.11(6X10) 


0.6012(88) 


0.981(14) 


0.9074(30) 


1.119(18) 


1.4074(58) 


-0.0711(21) 


(e) 


4.45(9)(20) 


0.808(33) 


1.1429(47) 


1.0004(17) 


1.311(11) 


1.4086(54) 


-0.1208(46) 



quark propagator which in turn produce disconnected diagrams for the hadronic correlators; 
denoting with ^U the quark propagator, one has Co,(d,«) (.xq, jo) = Lj(Tr{A^j75}Tr{A^575}) . 
This disconnected term is produced by 0{a) flavor symmetry breaking of the twist-term, similarly 



to that for the neutral pion in the light quark sector, cf. the contribution [ ]38[ ] and ref. [^] for a 
discussion of this subject. (The present case is however a bit different, since the disconnected term 
originates from direct mixing between c and s quarks, while in the neutral pion case it is due to a 
mismatch in the sign of the twisted mass for u and d quarks.) In this preliminary study we measure 
the Ds meson mass without taking into account the disconnected contribution. 

The simulation of the A^y = 2 + 1 + 1 theory is performed by a polynomial hybrid Monte Carlo 



algorithm (PHMC) []40|]. The structure of the algorithm goes along the lines indicated in []41|]. In 
table ^ we report the parameters of the first simulations performed. The parameters were chosen in 
such a way that the bare quark masses follow at full twist the proportions m„ : : : = 1 : 1 : 



5 : 60, see ( |7.3[ ) (of course renormalizations can alter these proportions). The point at K" = 0.1730 
is the one nearest to full twist where the pion mass is minimal: here ro/a ~ 3, nij^ ~ 400 MeV and 
the D meson mass is near to the physical value, but the kaon is heavier than in nature: roniK — 2.4. 
So the next step will be to decrease the s quark mass. Simulations on a 16^ • 32 lattice have been 
also started. 
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One interesting subject of investigation in perspective is of course the first order phase tran- 
sition at small quark masses in this new Nf = 2 + I + I setup, and in particular how it compares 
with previous studies with different gauge actions and/or quark flavor content (Nf = 2 and Wilson 
plaquette, DBW2 or tlSym gauge actions). The question arises, how the dynamics of the s quark 
can influence the phase transition and, more in general, the low-energy features of QCD. For the 
moment we can only state that no signs of metastabilities were found in the present runs, which 
could however be just too far from the metastable region. 



8. Conclusions 



We have investigated the phase structure of lattice QCD by probing different gauge actions, 
adding to the standard Wilson plaquette action a term proportional to the rectangular plaquette. 



see eq. ([2.2|). By varying the coupling ci which multiplies the rectangular plaquette term, one can 
interpolate between various actions and this allows to understand in more detail the properties of 
the phase structure, in particular how the strength of the transition depends on the additional term 
and how this influences the approach to the continuum limit. For the fermion action we always 
used Wilson twisted mass fermions at various values of the twisted mass parameter jj.. 

That even a small value of ci can already have a large impact on the phase structure is il- 



lustrated in figure [12| where we show the average plaquette value as a function of the hopping 
parameter K for three different actions, i.e. values of ci, namely ci = (Wilson), ci = — 1/12 (tl- 
Sym) and ci = —1.4088 (DBW2). As one moves K from the negative or positive side across the 



„ 0.56 - 



a 0.54 




G-O DBW2. |i=0.0 
□-0 tlSym. n=0.0 
O-O Wilson, n=0.013 




Figure 12: Hysteresis of the average plaquette value as K is moved across the critical point, for Wilson, 
tlSym and DBW2 gauge action at a 0.17 fm. 

critical point, where the PCAC quark mass vanishes, a hysteresis in the average plaquette value 
develops of which the size and the width is an indicator of the strength of the phase transition. We 
observe that both the width and the size of the gap in the plaquette value decreases considerably 
as we switch on c\ to ci = — 1/12 (tlSym action). Decreasing ci further down to ci = —1.4088 
(DBW2 action) still seems to reduce the size of the gap, but the effect is surprisingly small despite 
the large change in ci. Note that the results in figure 12 ai^e for a lattice spacing a ~ 0.17 fm that is 
roughly consistent for all three actions. 
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However, the results for the Wilson plaquette gauge action are at non-zero twisted mass = 
0.013, while they ai^e at zero twisted mass for the tlSym and the DBW2 action. Since the strength 
of the phase transition is expected to be reduced as one switches to a non-zero twisted mass, a true 
comparison at n =0 would disfavor the Wilson plaquette gauge action even more. Note also that 
the phase transition is a generic feature of dynamical Wilson fermion simulations independent of 
whether a twisted mass is switched on, i.e. the transition also occurs for standard Wilson fermion 
simulations. Another feature of the first order phase transition is that its strength weakens rapidly 
when the lattice spacing is made finer. This is illustrated in fig. ^for the case of the Wilson plaquette 
gauge action. 

A satisfactory setup for dynamical simulations with, say, Nf = 2 flavors of quarks would be 
to reach pion masses of about 250 — 300 MeV and a box size of L > 2 fm. At the same time, one 
should stay at full twist to realize 0(a)-improvement. From our preliminary results we find that 
for the tlSym action this can be achieved with a reasonable computer time at j8 = 3.9 on L/a = 20 
lattices. For smaller values of j8 we find that at pion masses of about 400 MeV large fluctuations 
appear, although no clear signs of metastabilities are visible. Although for the DBW2 action the 
situation might be somewhat better, the advantages of the tlSym action such as good convergence 
of perturbation theory and small scaling violations as found quenched, lead us to decide on the 
tlSym gauge action as the action of choice. 

For our present results from simulations using the Wilson plaquette and the DBW2 gauge 
actions we found that ;^PT describes the data rather well and allows to extract a number of low- 
energy constants of the effective chiral Lagrangian. We have also seen that scaling violations of 
both actions ai^e not very large when we compared physical results at different values of the lattice 
spacing. 

Very good news is that simulations with Nf = 2 + 1 + 1 flavors of dynamical quarks are per- 
fectly feasible with the twisted mass formulation of lattice QCD. Our first and preliminary results 
are very promising. The (PHMC) algorithm works nicely and, even more important, it seems that 
the tuning in this physically more realistic setup is quite manageable. Taken this fact together 
with the new developments of fast versions of the HMC algorithm to simulate dynamical quaite, it 
seems natural to start now realistic simulations with twisted mass fermions. 
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